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Abstract 



■ Transfer and decay dynamics of an exciton coupled to a polarization vibra- 

tion in a dimer is investigated in a mixed quantum-classical picture with the 
exciton decay incorporated by a sink site. Using a separation of time scales, 
it is possible to explain analytically the most important characteristics of the 
model. 

If the vibronic subsystem is fast, these are the enhancement of nonlin- 
ear self trapping due to the sink and the slowing down of the exciton decay 
for large coupling or sink strength. Numerical results obtained recently for 
the DST approximation to the model are quantitatively explained and new 



dynamic effects beyond this approximation are found. 

If the vibronic subsystem is slow, the behavior of the system follows closely 
the predictions of the adiabatic approximation. In this regime, the exciton 
decay crucially depends on the initial conditions of the vibronic subsystem. 

In the transition regime between adiabatic and DST approximation, com- 
plex dynamics is observed by numerical computation. We discuss the cor- 
respondence to the chaotic behavior of the excitonic-vibronic coupled dimer 
without trap. 
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I. INTRODUCTION 



The purpose of this paper is to study the interplay between a coherent transfer regime 
of an exciton and two processes leading to the loss of the linear character of the exciton 
transfer, namely trapping of the exciton at a sink site with a prescribed sink rate F and the 
coupling to intramolecular polarization vibrations. 

A lot of work has been done on exciton transfer theories during the last decades. Begin- 
ning with the microscopic treatment by Haken and Reineker and Grover and Silbey ^ 
a number of theories such as the Continuous Time Random Walk (CTRW) [Q, the Pauli 
Master equation (PME) @], the Generalized Master equation (GME) 0, the Stochastic Li- 
ouville equation (SLE) and the Haken-Strobl-Reineker (HSR) model (see and references 
therein) were developed and mainly directed to obtaining equations which describe the cou- 
pled coherent and incoherent motion of the excitation. The investigation of exciton transfer 
on relatively small molecular aggregates such as dimers or triades is of much interest for 
clarifying the applicability of exciton transfer theories to experimental situations such as 
described in 

Trapping of quasiparticles due to a sink site constitutes an important phenomenon in 
many molecular systems. In photosynthesis, for instance, an exciton in a harvesting an- 
tenna transfers its energy to a reaction center, where it can be trapped. Electron transfer 
processes then follow. Pearlstein and his coworkers were the first who recognized that the 
consequences of a sink on the energy transfer processes are different in the coherent and 
incoherent regimes (and references therein). Capek and Szocs |lO| pointed out the neces- 
sity of a transformation of the memory functions in the presence of a sink. They also gave a 
prescription for a proper inclusion of the sink into the HSR model. This found application 
e. g. in computer simulations of the excitation transfer in photosynthetic systems |lTT| , p!2| . 
Recently, a new form of purely coherent memory functions in presence of a sink was derived 
and shown to have important consequences for the excitation transfer |jl3[. New memory 
functions were also used to obtain characteristics entering the CTRW description in the 
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presence of a sink [ITl . 



The coupling between electronic and vibronic degrees of freedom in molecular and con- 
densed media is another basic mechanism influencing transfer properties of electronic ex- 
citations in these systems. The investigation of its consequences started from the polaron 
problem in solid states (see e. g. []T3[ and references therein; for exciton-phonon interaction 
see [|l6l) and continued with the study of the influence of the vibronic bath variables on the 
excitation transfer properties in the framework of the generalized Master equation P] and 
stochastic Liouville equation approaches 0. 

With the development of the theory of dynamical systems it has become attractive to 
analyze the implications of electronic-vibronic couplings employing concepts and methods 
of this field. Using such a dynamic system approach we study the detailed picture of the 
time evolution of a small number of relevant variables of the system, which are assumed 
to interact weakly with the environment. Recent experimental developments in the field of 
ultrashort time resolved spectroscopy (see e. g. ||T^) seem to make a direct observation of 
this time evolution possible in the near future. 

In the present paper we focus our interest on the exciton dynamics in a molecular dimer 
coupled to intramolecular vibrations. A remarkable feature of this model is the possibility 
of self trapping, i. e. unequal time averaged occupation probabilities on the two sites of the 
configuration. The easiest way to obtain this effect from a coupling to vibrational degrees 
of freedom leads to the two-site discrete self trapping (DST) equation [p!^-pO| which is a 
nonlinear but self contained equation of motion for the excitonic site occupation amplitudes. 
More sophisticated approaches take the dynamics of the vibrations explicitly into account 
by using a mixed quantum-classical description pT| , |22[| or by treating the coupled system 
quantum mechanically Another interesting point is the effect of dissipation on self 

trapping 

Although the influences of vibronic coupling and trapping on transport properties have 
been investigated separately in great detail, the combination of both, which can be important 
in the application of transfer theory, has rarely been addressed in the past. In a recent paper 



27| we have studied the interplay between vibrational coupling and trapping due to a sink 
site for a dimer and a trimer in the framework of the DST approximation. An extension of 
the results seems to be possible into two directions. One can either try to describe larger 
molecular aggregates within the same framework or one can try to improve the dynamic 
model which is used. In this paper we will adopt the latter way and take account of the 
dynamics of a particular vibronic variable, too. For this purpose the coupled excitonic- 
vibronic system will be treated in a mixed quantum-classical description |^ , ^ , which is 
justified whenever the quantum fluctuations in the vibronic subsystem are negligible. The 
model we use will be specified and developed further in section |I|, where we also indicate 
the modifications that lead to the DST equation. 

The explicit introduction of a vibronic degree of freedom makes not only the dynamics of 
the system more complex, it also increases the number of initial conditions that have to be 
specified in order to uniquely determine a solution. Motivated by the experimental situation 
e. g. in photosynthetic units, the creation of the exciton will always be assumed at the site 
without sink. For the vibronic degrees of freedom we will consider various possibilities. As 
we shall see, the initial conditions can crucially influence the dynamics and it is not easy 
to overlook all the dynamical regimes the model is capable of from numerical computations 



only. Therefore we devote section |Tl| to an analytical treatment of the system assuming some 



separation between the different time scales. In particular, we will generalize the analysis 



of the fixed points for the excitonic-vibronic coupled dimer without sink p2[ to the present 
system. This will provide a global picture of the phase space which we shall further support 
in section |^ by discussing the time dependence of the total occupation probability and the 
relative site occupation probabilities for some specific solutions. 

The results of section ^TT| will also allow for a quantitative understanding of the DST 
dimer with sink. One of our aims is to study what is inherited in our more complex model 
from this approximation and which dynamic effects are beyond it. We will therefore briefly 
recall the relevant numerical results from |^ and discuss them in the light of our findings 



in section III A 



Since the mixed quantum-classical description we are using in this paper can be justified 
best for small oscillator frequencies, we pay much attention to this adiabatic regime, too. 
As we shall see, the assumption of slow vibrations is just the antipode of the DST case and 
therefore new dynamic effects different from DST results can be expected from it. 



II. DESCRIPTION OF THE MODEL 



A. The Hamiltonian 



We consider the dynamics of an exciton moving on a molecular dimer. At each of the 
two monomers the exciton is allowed to interact with an intramolecular vibronic degree of 
freedom. Thus, the Hamiltonian of our model contains excitonic, vibronic and interaction 
parts denoted by ifcxc, H^ih and i^int, respectively: 

H = Hcxc + H^ih + -f^int, (1) 

Hcxc describes a two site model 



with?!, m = 1,2. c„ is the probability amplitude of the exciton to occupy the n-th molecule 
and Vnm the transfer matrix element due to dipole-dipole interaction. In a standard two 
site model, the are real quantities and correspond to the local site energies of the exciton. 
Here, we allow €2 to contain a negative imaginary part in order to describe the decay of the 
exciton on the sink site 2. Since we are not interested in the effect of a site energy difference 
on the exciton dynamics here, we set 

r 

ei = 62 = (3) 



which is equivalent to the extended sink model for the exciton decay introduced in [10 



on 



the density matrix level. This model has been shown to solve the problem of a consistent 
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description of exciton trapping at a sink meeting basic physical requirements such as positive 
occupation probabihties. 

The vibrational part i^vib is taken as the sum of the energies corresponding to intramolec- 
ular vibrations at each of the monomers for which we use the harmonic approximation 

^vib = E^(Pn + ^X)- (4) 

Here qn-.Vn and ujn are the coordinate, the canonic conjugate momentum and the frequency 
of the intramolecular vibration of the ra-th molecule, respectively. 

The interaction Hamiltonian takes into account that the exciton energy depends on the 
molecular configuration of the monomers which is expressed by the coordinates g^. Using a 
first order expansion in g„ one has 

-f^int = XI 7ngnC* Cn, (5) 
n 

where 7„ are some coupling constants. In order to restrict the number of free parameters as 
much as possible we will assume that the dimer is symmetric except for the additional sink 
term on site 2, i. e. c^i = = i^, 7i = 72 = 7 and V\2 = V21 = —V. 

In what follows we will use a mixed quantum-classical description of the dynamics, 
i. e. we treat the vibronic degrees of freedom in the classical approximation while retaining 
the quantum wave function for the excitonic two site system. This approximation can be 
justified over a finite time range which increases as the oscillator frequencies and coupling 
constants decrease. When it exceeds the life time of the exciton the mixed quantum-classical 
picture describes correctly the decay of the excitation. 

Using units with ^ = 1 we obtain from ([l|)-(|^) the equations of motion 



id/dtci = 7giCi — V C2 (6) 

i d/dt C2 = (^-i^ + 7^2) C2-VC1 (7) 

d/ dt g„ = p„ (8) 

d/dtpn = -Uj'^qn-'y\Cn\'^- (9) 
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B. Reduced Equations of Motion 



For a numerical investigation the equations are well suited and we have integrated 

them in order to obtain the results that will be presented in section |I^. The analytical treat- 
ment of section |T|, however, requires to reduce the number of variables and free parameters 
as much as possible. Therefore we will now rewrite the equations of motion (^)-(y) using 
appropriate dimensionless variables and parameters. 

The excitonic subsystem can be described by a point on the Bloch sphere which is usually 
given in Cartesian coordinates. Here we prefer to parametrize the Bloch sphere in spherical 
coordinates i?, 6 and 0. They are defined using the density matrix of the two-site system 

Pmn = CmCn (n,m = 1,2) 

R ■= PU + P22 

cose R := P22 - pu {0<e < 7i) (10) 

e^'^sin^ i? := 2pi2 (-vr < < vr) . 

Due to the trapping of the exciton the radius of the Bloch sphere R{t) which is the total 
probability to find an exciton on either of the two sites is not constant but a monotonically 
decreasing function of time with R{0) = 1. 

Besides the total occupation probability R the difference of the occupation of the two 
sites is of interest. It is determined by the angle 6 since we have 

2 l-cos6' 2 l + cos6' 
|cir = ^ R |C2|' = R. (11) 

The phase has no direct physical interpretation. We note that is not well defined at 
the points ^ = and 6 = n. This can be circumvented by directly considering the time 
dependence of p at these points and will not affect the following. 

When deriving the equations of motion for the new variables from (^-(§) one observes 
that only the difference q2 — qi couples to the excitonic degrees of freedom. Therefore we 
can introduce a dimensionless difference coordinate and the conjugate momentum 

8 



Q := VV{q, - gi) P := -^{P2 - Pi) (12) 

and reduce the number of independent variables in this way by two. The reduced equations of 
motion for the remaining five variables are obtained after the introduction of a dimensionless 
time 



T = 2Vt (13) 



and dimensionless parameters 



■= r-=— ■=— (U) 

describing the strength of the electronic-vibronic coupling, the frequency ratio of the two 

interacting subsystems and the strength of the sink, respectively. We find 

R = -g{cos9 + l)R (15) 

6 = g sinO + sin(f) (16) 

^ = cot 9 cos (p — ^J2p r Q (17) 

Q = P (18) 



P = -r^Q-^p/2rRcose. (19) 
In these equations (') denotes d/dr. 

C. The DST Approximation 

One standard way to simplify the dynamics of excitonic-vibronic coupled systems is 
to assume that the vibronic degrees of freedom instantaneously adapt to the state of the 
excitonic subsystem and always remain in the ground state prescribed by it. Applied to 
(0)"(S)' ^^^^ assumption results in the DST equations mentioned in the introduction. 

Within our effective dynamic model (|l5|)-(p!9D the DST approximation can be justified 
assuming a separation of time scales. The time scale for the (free) transfer of the excitation 
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between the two sites of the dimer has been normahzed to 1 when the equations of motion 



were written using the dimensionless time r (13). Another relevant characteristic time of 



our system is the period of the oscillator which is ~ 1/r. Now we assume that the vibronic 
degrees of freedom are much faster than the exciton, i. e. r ^ 1. In this case the oscillator 
coordinate completes many cycles during a time of the order 1 which is relevant for the 
slow excitonic subsystem. Consequently, dynamic self averaging over the fast oscillator 
coordinate occurs and Q in (|17|) can be replaced by its time average. For the harmonic 
oscillator (Q, P) this average is given for arbitrary amplitude by the oscillator ground state 
which is determined by the state of the excitonic subsystem and which represents at the 
same time the only fixed point of the oscillator dynamics, i. e. formally we can introduce 
the DST approximation by requiring quasistationarity in the vibronic variables Q = P = Q. 
We obtain for the mean value of Q 

Qbst = --\I?:Rcos9. (20) 
r V / 

and substitute it into ( |T7| ) which is replaced by 

(p = cot9 cos(f) — p R cos9 . (21) 

Together with (|T^) and (|T^ this equation governs the DST dynamics of our model. 

An analytical justification for the averaging procedures applied in this and the following 
sections can be given using mathematical tools that were developed in the theory of nonlinear 
differential equations (see e. g. p8[) and will not be discussed here. Instead we confirm the 
resulting equation (^) by the numerical simulations for r ^ 1 presented at the end of this 
paper. 

The derivation of the DST approximation using fast oscillator dynamics is questionable 
although it seems straightforward within the mixed quantum-classical description to which 
this paper is confined. However, the assumption r ^ 1 means that the mixed quantum- 
classical description itself looses its justification and should be replaced by a full quantum 
treatment. More consistent ways to address the validity of the DST limit are based on 
dissipation due to a quantum heat bath |^ and therefore beyond the scope of our model. 

10 



The DST equation is known to reproduce at least qualitatively some remarkable features 
which the full dynamic system displays for an arbitrary value of r, e. g. the bifurcation in 
the phase space for p = 1 ||18[, the resulting possibility of self trapped solutions for coupling 



strengths above this value |19] and the possibility of dynamical chaos when an external 



perturbation is applied |]2D|. We can therefore consider the results obtained in within 
the DST approximation for a dimer with sink as a guiding line for the effects that can be 
expected in the present more complete treatment. 

III. QUASISTATIONARY DECAY MODES 
A. Fixed Points for Quasistationary Total Occupation 

Beside the free exciton transfer time and the oscillator period a third relevant time 
scale controls the decay of the excitation. As we shall see it is not necessarily given by 
the inverse sink rate 1/g. In the present section we will assume that the exciton decay is 
the slowest process in the system and discuss the resulting quasistationary solutions of the 
equations of motion ([T5|)-(p!9D. In the special case of the DST limit these quasistationary 
solutions are most transparent, while the opposite case of an oscillator which is even slower 
than the excitation decay will be discussed in the next subsection and results in different 
quasistationary solutions. 

The further analysis is motivated by the observation that a self contained equation for 
the decay of the total occupation probability R{t) would be obtained if the function 9{R) 
was known. This is the case, in particular if the system remains in a fixed point of the 



equations ([T6|)-(|T9D while R varies slowly enough to be considered as an adiabatic parameter 
of these equations. Then, according to (0), either the sink rate g is very small or the 
quasistationary state has a site occupation difference which is strongly biased towards the 
site without sink cos^^ ~ — 1. 

In the following we discuss the fixed points of the equations (p!6D-(p!9D with R treated 
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as an adiabatic parameter. The equations ([T8|) and (|T9D yield for a fixed point P = 
and Q = Qbst, i- e. the location of the fixed points is the same for the system with the 



full oscillator dynamics included and for the DST equations. The equations (|1^) and (|T1 
allow for two different pairs of fixed points on the Bloch sphere classified in what follows as 
detrapped and self trapped states. The stability exponents for any of these points can be 
obtained from a linearization of the equations of motion (|T6|)-(p!9D around the fixed point 
which yields the characteristic equation 



cos^ 



= (A^ + r^) ([A + sin(/) cote][X - g cosO] + j - r^pR cos(f) sinO . (22) 



(A) Detrapped states 
For sufficiently small sink rate g < I we obtain two fixed points at 



sin (j) = —g, cos0 = ±yl — (23) 

cose = 0. (24) 

Because of ( p4|) the occupation probabilities for the two sites are the same and we call 
the fixed points detrapped states. The point A+ at cos(f) > can be considered as 
a generalization of the bonding state in the system without sink whereas the point at 
cos (p < corresponds to the antibonding state. The stability exponents of these fixed points 
are given by 



- 2 + COS^ , 

A^ = ± 



— cos^ I 



2 



+ r'^pR COS0 . (25) 



For A"*" the argument of the square root is always positive and the pair of stability exponents 
corresponding to the negative sign in ( P^D is purely imaginary. The other pair consists of 
two imaginary or two real exponents with opposite signs such that the fixed point A"*" is 
elliptic if {pRY + g'^ < I and unstable hyperbolic otherwise. 



12 



For A~ the stability is determined by the argument of the square root in (]25|). If it is 
positive, the point is a stable elliptic center and this is always the case in the adiabatic regime 
r — >■ 0, in the opposite DST case r ^ cxo or for arbitrary parameters at large times since 
then R 0. Only for p sufficiently large the argument of the square root may temporarily 
be negative. The stability exponents then acquire real parts with opposite signs and the 
point A~ renders unstable hyperbolic. 

For the DST case r ^ 1 one pair of stability exponents which is given by A = ±ir 
corresponds to the fast oscillations around the DST solution whereas the other pair of 
exponents 



describes the stability of the DST solution itself. 

We note, that the positions of the two fixed points do not depend on R and are 
thus constant in time. For the full system with R time dependent they represent therefore 
special time dependent states in which the distribution of the excitation over the two sites 
is constant, just the total occupation decreases exponentially at a rate g 



When these fixed points are stable, a state prepared in their vicinity will remain there and 
decay at a mean rate g with some oscillations superimposed. 



For a strong sink or strong nonlinearity (pR)"^ H-fi'^ > 1 there exist two other fixed points 
with biased site occupation probabilities, i. e. self trapped states: 



= cos (f){pR — cos 0) 



(26) 



R=-g R. 



(27) 



(B) Self trapped states 



sin 6 



1 




(28) 




sin (j) = —g sin 6 



cos (f) = pR sin 6 



(29) 
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The existence of these two fixed points corresponds exactly to the range of parameters for 
which the fixed point A"*" is hyperbohc and for (pR)'^ + g"^ = 1 the points A+ and merge 
into a single one. So we have established a generalization of the pitchfork bifurcation of the 
system without sink |^ . The difference is that the bifurcation parameter no longer depends 
exclusively on the coupling strength p. Rather it contains along with p the strength of the 
sink g and the total occupation probability R which is a function of time. U g < 1 the 
fixed points will disappear for large times when R ^ and we can speak of a dynamic 
bifurcation. On the other hand, when g > 1 the detrapped states A do not exist anymore 
and there is no bifurcation in the course of time. 

Another crucial difference to the sink less dimer has to do with the character of the points 
B. Without sink they are stable elliptic centers and in their vicinity there exist solutions 
which remain self trapped for all times. In the present case the stability exponents have to 
be determined from the quartic equation 

= (r2 + A2)([A cos^]^+p2i?^) - {rpR sinOf , (30) 

which does not allow for an easy solution. We shall see later on, that the fixed points B^ due 
to their time dependence do not represent quasistationary solutions unless the bifurcation 
parameter sin^^ ^ = {jpRY + g^ is far enough above the bifurcation value 1. Therefore we 
simplify (|30|) under the assumption sin^ ^ 1 and drop the second term. Then, there are 
two pairs of solutions for the stability exponents. One of them is purely imaginary A = ±ir 
and the other one contains a real part as well 

\ = gcos9±ipR. (31) 

Again, for r — > oo the first pair describes the oscillations around the DST solution and 
the other one the stability of the DST solution which can alternatively be obtained from a 
linearization of (|T6|) and (|21|) without further approximations as 

A = 5'cos6' ± i I cos6'|pi?. (32) 
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Note that ± in ( ]28|) stands for the two different fixed points, whereas the ± in ([3l| ) and (|32D 
corresponds to two different stabihty exponents A of the same fixed point. 

The character of the fixed points is determined by the real part of A. For with 
an occupation bias towards the sink site (cos^^ > 0) we have a repeller whereas the point 

with a low occupation probability at the sink site (cos^ < 0) is a stable attractor. 

It is less clear than for the detrapped states that the fixed points B which were obtained 
under the assumption of a constant total occupation have some interpretation for the full 
system since their location does depend on R{t). However, we will show now that at least 
B^ can represent an attractor and a quasistationary solution for the complete equations of 
motion (p!5D-(p!9|) when we are sufficiently far away from its threshold of existence, i. e. when 
the bifurcation parameter is sufficiently large 

sm^g= ] «1. (33) 

For this purpose we have to show that the change in the position of B~ is much slower than 
the relaxation towards this fixed point. The latter occurs on a time scale given by (|3lD as 
1/g, while the oscillations around the fixed point have a period 2n/r and will be averaged 
out provided this time is small enough. 

The velocity of the fixed point location can be estimated after inserting ([T5|) into (|28D 
and (|2|). We obtain 

2sin^6'(l + cos6') 



cos 6* I 



sin^ 9 



\e\=g{pRy 
g{pRf . 

c 

< - sin^ e 
2 

ill / „,2sin^^(l + cosi 
101 = {gpRY 



COS( 

sin^^ 



g^pR^: 4 



2 

2 

and conclude that the position of B~ changes always slowly when sin^ 9 1 and in par- 
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ticular the relaxation towards the attractor is much faster. If moreover g/2 sin^ ^ r the 
oscillations around the fixed point average out in ( [T5| ) and we can consider (-R(t), ^^[i?(r)], 
0[i?(r)]) as an attractor for the full system. This condition is automatically satisfied when 
r > 1 (e. g. in the DST regime), but in the adiabatic regime r <^ 1 the assumption of quasis- 
tationarity of B~ may not be valid or restricted to a short interval in time. In this case the 
following does not apply and has to be replaced by the discussion in the next subsection. 

Once the system is close to the attractor B~ the excitation decay is approximately 
governed by the equation 

R=- , ^J?' , R. (34) 



If the sink rate is dominant pR < g we obtain from (|3^) the interesting effect that the decay 
rate decreases for increasing sink rate g. The reason is that the location of B~ is strongly 
shifted to the site without sink such that the probability to find the exciton on the sink 
site becomes very small. This behavior has previously been studied without coupling to 
vibrations and was termed "fear of death" effect in some publications ||rT|. When we have 
pR/g which is always the case for r — oo (^) implies an exponential decay of the 
exciton with the rate l/2g. From ( [27| ) and ( |5^ we conclude that the fastest quasistationary 
decay of the excitation is realized for small electronic- vibronic coupling and a sink rate g = 1- 

Self trapping on the site with sink is not immediately destroyed when the sink becomes 
effective. Though the state is a repeller, the solution can oscillate with slowly increasing 
amplitude around this fixed point provided that g is not too large. In this case B+ can control 
the dynamics for some finite time which then leads to an enhanced excitation decay. We do 
not want to pursue this possibility further, since we restrict our discussion to an excitation 
initially localized outside the sink and in this case the point B+ is always irrelevant. 

To end this section we would like to recall some of the numerical results obtained in ||2^ 



for the DST dimer. There we found that the system relaxes to a quasistationary self trapped 
state provided that the nonlinearity was sufficiently large. Moreover, it was demonstrated 
that the sink supported the tendency to self trapping and that the self trapped state disap- 
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peared for large time when the sink was not too strong. For strong sink we showed that the 
hfe time of the exciton was the larger the stronger the sink was chosen. Now we are able 
to interpret these qualitative observations as the dynamical manifestations of the point B~ 
and can give them a quantitative formulation using the results of this subsection. 



B. The Adiabatic Regime 

In the previous subsection we considered quasistationary dynamics of the excitation 
under the assumption that the decay of the exciton represents the slowest process in the 
system. This assumption breaks down when the system is in the adiabatic regime r ^ 1. 
Then the decay is an essentially nonstationary process leading to the disappearance of the 
exciton before the oscillator has gone through a large number of cycles. 

Opposite to the derivation of the DST equation we can now assume that the exciton 
completes many oscillations on the time scale of the oscillator. Again we want to exploit 
this fact by averaging over the variables of the fast subsystem and replacing them by their 
mean value, but in contrast to the harmonic oscillator of the DST case, the equations of 
motion for the excitonic variables have two fixed points for a given Q which are obtained by 
setting the 1. h. s. of eqs. ([T^) and (|16D to zero. This results in 



sin =—(7 sin ^ cos = y 2pr Q tan6' (35) 

and combining these two equations we find 

c/2 sin^ e + 2pr'^ tan^ = 1. (36) 

The latter equation determines the location of the fixed points which depends parametrically 
on Q. The stability exponents can be given in the form 



\ = gcose±i^2pr\Q/ cose\. (37) 

Eq. ( P^D is a biquadratic equation in cos 6*, i. e. the two fixed points have opposite signs for 
cos^^. Due to (|37D this means that one of them is a stable attractor, the other one a repeller. 

17 



Instead of writing down the explicit solution of (^) which is quite a lengthy expression 
though easily found we would like to mention two limiting cases. 

First we note, that as g —>■ the two fixed points approach the well known lower and 
upper adiabatic states of the system without trap (see e. g. [p3|) which are given by 



cos e = ±—^L^— cos </) = ±1 . (38) 

The second important limiting case corresponds to a strong localization of the exciton on 
one of the two dimer sites (cos 6* — ±1). Assuming g'^ + 2pr'^ ^ 1 the solution approaches 

In particular, for a very large sink rate g the dependence on Q disappears and we have 

sin 6* = - < 1 sin = -1 . (40) 

In order to be able to replace the time average for the excitonic variables by the discussed 
fixed points we have to assume that the time scale given by ( P7| ) for the relaxation is 
sufficiently short, or that the oscillations around the fixed point do average out. Due to 
the nonlinear equations for the excitonic variables and in contrast to the derivation of the 
DST equation we have in this latter case to assume that the amplitude of the oscillations is 
small, i. e. the system has to be prepared close to one of the adiabatic states. If so, even the 
adiabatic state with cos^ < 0, which is in fact a repeller, can be considered quasistationary 
for some limited time. 

A self contained equation for the decay of the total excitation probability R{t) can be 
derived under the assumption that the exciton is located close to the attractive quasistation- 
ary state. Then cos^ may be replaced by the value prescribed by the oscillator coordinate 



Q according to (pq ) and is constant provided the oscillator dynamics is sufficiently slow to 
be completely disregarded during the life time of the exciton or if the position of the fixed 
point does according to (^) not depend on time due to a strong trap. Either case leads to 
an exponential decay of the excitation. 



IV. THE TIME EVOLUTION OF THE SYSTEM 



A. Parameter Regions and Initial Conditions 

In the previous section we have discussed possible quasistationary states of our system 
without indicating, if and when these states will be reached by a particular solution of the 
equations of motion. Now we turn to the investigation of the time evolution of particular 
solutions starting from specified initial conditions. 

As mentioned in the introduction we will always consider an exciton which is created at 
the site without trap, i. e. we set initially 

ci(0) = l C2(0) = 0. (41) 

The initial conditions for the vibronic degrees of freedom which we have chosen are meant 
to take into account different physical possibilities to prepare the excitation and to provide 
enough variety to estimate the degree to which the excitonic variables depend on the details 
of the oscillator initial state. The resulting solutions will be referred to in the following way: 

1. Bare exciton: The first initial state we consider corresponds to a creation of the exciton 
on the first molecule by a very short light pulse. During the optical excitation from the 
ground state there is not enough time for the local vibrations to accommodate to the 
creation of the exciton. This means we assume the vibrations initially in their ground 
state without exciton q'j(O) = 0, Pi(0) = 0. The total energy for this initial condition is 
0. 

2. Polaron: The second possibility is to assume a slow excitation such that initially the 
vibrational degrees of freedom arc already relaxed to their new ground state with exci- 
ton. This is the initial condition which would be implied by the DST approximation: 
gi(0) = —^ju?, ^2(0) = 0, pj(0) = 0. The total energy of the polaron is —7^/20;^, i. e. 
lower than the energy of the bare exciton. 
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3. Polaron with additional vibrational energy [Polaron (-) and Polaron (+)\. The differ- 
ent initial energies make a direct comparison between the bare exciton and the polaron 
difficult. Therefore we have taken into account a third possibility for the initial condi- 
tion which is not directly related to a particular way of preparing the exciton. Again 
we chose the configuration coordinate of the vibrations in the minimum of the poten- 
tial after the exciton has been created. We supply, however, an initial momentum such 
that the total energy is as for the bare exciton case. For this momentum we have 
two different possible directions i. e. the polaron (±) is specified by gi(0) = — 7/0;^, 
g2 = 0, pi(0) = ±7/c^, P2(0) = 0. 



Before we explore the different types of solutions using the results of section |T| as well 
as numerical simulations we would like to briefly mention some limiting cases for which 
analytical solutions to our problem can be given. For a sinkless DST dimer [g = 0) exact 



solutions in terms of elliptic functions were found in |T9|. Starting from a state completely 
localized on one of the dimer sites they turned out to be self trapped for p > 2. Recently, for 
certain parameter regions and g = exact solutions with the oscillator dynamics explicitly 
taken into account in a mixed quantum-classical description were also obtained pT| . 

For the linear dimer {p = 0) with a sink the equations of motion can be integrated and 
yield for the total occupation probability 

R(r) = (l-G' + ^!^e+v^^ + 9l±^e-^A e"^- (42) 



with G := — . For the full set of equations of motion (|T5D-(|T9]) we have analytical 



solutions for exceptional cases only such as the quasistationary states of section Hi A . 
Even for the DST approximation with sink no exact solutions are known in general. 

We have performed a numerical integration of the coupled system of equations (^-(^ 
for the different described initial conditions. We display results for the total occupation 
probability R{t) and the relative site occupation difference expressed by cos 6'(r) for various 
values of the oscillator frequency ranging from the high frequency (DST) limit in figs. to 
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the deeply adiabatic region in fig. We restrict the sink rate and the vibrational coupling 
to three representative cases: 

• weak sink = 0.1 / weak coupling p = 1: fig. |1| 

• weak sink g = 0.1 / strong coupling p = 3: figs. ^, ^, |^ and ^ 

• strong sink g = 3 / strong coupling p = 3: figs. § and ^ 

For each set of parameters the results for the different initial conditions will be displayed in 
the same graph. They can be distinguished by the different line shapes annotated e. g. in 
fig. 1(a). 



B. Time Evolution in the DST Approximation 

Using the results of section |111 A| we can obtain a quite satisfactory description of the 



time evolution in DST approximation that agrees with our numerical findings reported in 
| P7| . We have to distinguish three different cases with respect to the parameters p and g: 



(I) Nearly Linear Regime + < 1 
In this case throughout the whole time evolution the only fixed points present are the 



stable elliptic centers A from section [III A| . The relative site occupation difference cos 9 will 
therefore oscillate with a mean value cos ^ = 0. The decay of the total occupation probability 
is then approximately given by (pT]): 

R{t) ~ , (43) 



and is therefore very much like in the case of the linear dimer 

An illustration for the described behaviour is provided by fig. |1|, which is with p = 1 
and g = 0.1 at the fringe of region (I). Since the transition between the parameter regions 
is smooth we find in fig. |l|(a) a straight line indicating an exponential decay with some 
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oscillations superimposed. The mean decay rate obtained from the figure is in good cor- 
respondence to ( ^31) very close to 0.1 and the period of the oscillations in fig. |l](b) is very 
close to 27T. This is the value for the free transfer of the excitation and corresponds to 
the asymptotic value for the stability exponent of the points A±. However, the solution is 
actually not in the vicinity of one of these points. Rather it oscillates with a large amplitude 
and can therefore not be expected to be correctly described by a linearization around a fixed 
point. For instance, the time dependence of the stability exponent ( p6|) is not refiected in 
the solution. 

(11) Weak Sink g < 1 and Strong Coupling g"^ + > \ 

In this case there exists the attractive fixed point B~ from section 3.1 when the system 
starts its evolution at R{t = 0) = 1. As a numerical example consider the DST curve of fig. 
which is the thick gray line. The system approaches the attractor after a time r ~ 1/(7 and 



then decays on it according to (p^). This is in general a nonexponential decay which is very 
much different from p3|). If we assume strong nonlinearity p ^ 1 (0) can be approximated 
and results in 



R{r) = ^l-f,r. (44) 

There will be oscillations around this mean behavior with an amplitude decreasing as the 
attractor is approached. The frequency Vt of these oscillations is given by the imaginary 
part of the stability exponent in (0) and decreases approximately as ~ pR for strong 
nonlinearity. Indeed, the oscillations around the mean in fig. ^(b) have a period which can 
be seen to increase starting from T ~ 2.4 which is close to the value 2.2 obtained from 
the imaginary part of (^). Then the oscillations die out at r ~ 30 thus confirming the 
attractive character of the fixed point B~. 

When the total occupation has decreased such that g'^ + (pR)'^ ~ 1, the attractor B~ 
does not exist anymore and the system will start oscillating with equal mean site occupation 
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probabilities around one of the fixed points as in (I). The time tq for the crossover from 
the algebraic decay (^4|) to an exponential behaviour with decay rate g is approximately 
given by 

n2 - 1 

ro = ^ . (45) 

9 

Of course the crossover does not occur instantaneously, rather there exists a time interval 
close to the tq where we do not have a good description for the system. We recall that the 
assumption of a quasistationary decay on the point B~ required the system to be sufficiently 
far away from the dynamical bifurcation (|33|). This means that the crossover will actually 
start a little earlier than predicted by (|45|) as can be seen in fig. ^ where the crossover 
according to (^5]) should be at r = 80 but occurs actually at r ~ 70. It means also that the 
dynamics for + ~ 1 will be hardly different from the case (I), where an example for 
this transition behavior was already described with fig. ^ 

( III) Strong Sink g > I 

Here the attractor does exist throughout the evolution of the system. If the nonlin- 
earity is very large there might be initially a nonexponential behavior as in (II), but this will 
turn into an exponential decay as soon as the total occupation has decreased sufficiently for 
g ^ pR. Then one has from 



R{t) ~ e-(^-v^)" . (46) 

Under the assumption p ^ g which is, however, not satisfied in the numerical example fig. 0, 
the approximate crossover time is obtained from (|4^) as 



To = — - — • (47) 

The DST solution represented by the thick gray line in fig. ^ relaxes after a very short time 
to the attractor whose initial and final position is marked by dotted horizontal lines. 
Due to the larger g compared to fig. |^ the initial oscillations around B~ can hardly be 
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observed. The crossover to constant relative site occupation probabilities and exponential 
decay occurs at r ~ 20. 

We note that in all the three cases we discussed the limiting behavior for large time r 
and small total occupation R agrees with that of the linear dimer as given by (^2]), i. e. it 
is exponential. This is natural, since the second term in the DST equation (^Tj) which due 
to the excitonic-vibronic coupling is then negligible, but it is in contrast to the full system 
of coupled excitonic-vibronic equations ([T5D -(|l9|) from which it is obvious that the oscillator 
keeps influencing the dynamics of the exciton for all times. Just the feedback from the 
exciton disappears as it decays. 



C. Deviations from the DST approximation at large but finite oscillator frequency 

The three different scenarios for the time evolution in DST approximation which were 
described in the previous section remain valid for the full system with the oscillator frequency 
not too low, since the fixed points of p.llA| still represent quasistationary decay modes. Using 



the numerical results displayed in fig. we will demonstrate this, discussing at the same 
time deviations from the DST solutions. 

The degree of deviations from the DST solution will depend on the value of the parameter 
r and on the initial conditions for the oscillator. In particular for intermediate oscillator 
frequency it can be expected that the DST approximation describes the actual solution the 
better, the closer to it the initial condition for the oscillator is chosen. Indeed, the polaron 
(0) (dashed thick gray line), which is prepared in a DST state, cannot be distinguished at all 
from the DST curve in the plots for high oscillator frequency r = 10 (fig. and follows it 
very closely in the plots for the intermediate frequency r = 1 (fig. ^). For small nonlinearity 
p = 1 (fig. |1|) the same holds true for the other three solutions which are prepared with higher 
energy than the DST solution and in this case the initial conditions have no crucial influence 
on the dynamics down to the intermediate oscillator frequency r = 1 (not displayed). 

A systematic, though small deviation from the DST solution can be observed for stronger 
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nonlinearity in the figs. 0, ^and ||. Here, a self trapped state is approaclied by tlie polaron(±) 
and tlie bare exciton solution as described for the DST case, but beside the familiar slowly 
decaying oscillations we find also oscillations of higher frequency which do not disappear 
completely. This behavior was to be expected from the stability analysis of the fixed point 

for the complete system, where we found from ( |50D beside the stability exponents of the 
DST solution a pair A = ±ir describing fast oscillations. However, the location of the fixed 
point and the center of the oscillations in figs. |^ and |^ are not exactly the same. The full 
dynamic model tends to oscillate around a state which is even more localized than predicted 
and consequently it decays slightly slower. Moreover, the period of the fast oscillations is of 
the order of 27r/r but does not quite agree with this value and is actually close to half of it. 

The reason for these deviations can be seen in the fact that a linearization of the flow 
around the fixed point is justified for small amplitudes only - a condition not matched by 
the vibronic variables, since the oscillator was prepared in a state of high energy. Moreover, 
after r ~ 20 when the exciton has almost ceased to exist a description using fixed points of 
the coupled excitonic- vibronic system - though it still provides a fairly good picture - seems 
to be counter intuitive since we have seen that the interaction in this case is one way only: 
the freely moving oscillator represents an external perturbation to the excitonic subsystem. 
We will come back to this point when we discuss the adiabatic case in the next section. For 
the time being it is sufficient to note that the deviations shrink for growing r as it can be 
seen by comparing figs. ^ and |^ and have completely disappeared on the scale of the plots 
for r = 100 (not displayed). 

Unlike in the figures discussed so far, a qualitative difference in the behavior of the three 
solutions prepared with total energy is observed for p = 3 and g = 0.1 when they are 
compared to the low energy DST and polaron (0) solutions. In fig. ^(b) beside the latter 
two the bare exciton is shown for which the oscillations around the initially existing self 
trapped state B~ are so large that this state is hardly recognized at all. It disappears at 
r ~ 10 and thus much earlier than for the DST and polaron (0) case. The other two polaron 
solutions which are not displayed resemble the bare exciton. 
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This behavior can be understood from the fact that the system with g = 0.1 is very 
similar to the sink less case which has been shown to be strongly chaotic for r = 1 and p 
above the bifurcation value 1 |2^. Due to the chaos, the system explores the energetically 
accessible phase space very fast and this is reflected in the strong and irregular oscillations 
of the relative site occupation leading to a rapid exciton decay. 



D. Time Evolution in the Adiabatic Case 



The conclusion of section |III B| was that there is an exponential decay once the system 
is close to the attractive adiabatic state. However, there are two important limitations to 
this conclusion. 

First, the initial state of the exciton has to be close to one of the adiabatic states. We 
consider the initially completely localized state of the exciton cos^^ = —1, i. e. the initial 
oscillator coordinate should correspond to a strongly localized adiabatic state. According to 
(^) this is the case if + 2pr^Q(0)^ ^ 1. For the bare exciton Q{0) = this is the case 
for large sink rate g and then the total occupation will decay close to the linear dimer (^21) 
independent on the actual strength p of the coupling. For the polaron we have 

/ + 2pr2g(0)2 = /+p2 > 1 (48) 

and the exciton can be close to an adiabatic state even for small sink rate provided the 
coupling to the oscillator is strong enough. The total probability decays in this case as 

R{t) ~ e"^" (49) 

which is initially very close to what is predicted from the quasistationary decay mode B~ 
(p^). The largest possible decay rate is according to ( ^Qf ) l/4p and it is realized for p = g. 

The second condition for an exponential decay of the occupation probability is that the 
oscillator dynamics is actually sufficiently slow to be completely disregarded during the life 
time of the exciton. According to (EUI) this means 
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but the restriction of the exciton to one of the adiabatic states prescribed by the oscillator 
is justified whenever r ^ 1, and this can be a much weaker condition. So if in the adiabatic 
regime the condition (^) is not satisfied, no self contained equation for the decay of the 
exciton is available. 

This is the situation in fig. ^ The parameter r = 0.1 is sufficiently small for the applica- 
tion of the adiabatic approximation. Consequently, in part (b) of the figure the polaron (-) 
solution can be seen to follow the evolution of one of the adiabatic states, namely the ener- 
getically lower state obtained from the solution of (|36D . Initially, some decaying oscillations 
around the adiabatic state can be observed which are in good agreement with the stability 
exponents (^71). When the adiabatic state enters the region cos 6' > it becomes a repeller 
and one observes increasing oscillations around it until the variable Q has completed one 
full period at r ~ 80 and the relaxation to the attractor starts again. 

The same behavior can be observed for the other two polaronic solutions. The excitation 
decays rapidly as soon as the exciton is driven by the oscillator to the sink site. Therefore 
in this case the life time of the excitation is basically determined by the frequency of the 
oscillator and its initial conditions. Since the polaron (-) has an initial momentum which 
is directed towards increasing polarization, the exciton remains for a long initial period 
localized on the site without sink. This period is shorter for the polaron (+) which has a 
momentum towards decreasing polarization and consequently the polaron (-) has a longer life 
time. The polaron (0) has no momentum at r = and decays initially at a rate in between 
the other two polarons. Due to the lacking vibrational energy the oscillator coordinate for 
the polaron (0) changes its position only very slowly such that the polaron (0) is the longest 
living solution. 

Comparing in fig. |^(a) the polaron (0) to the exciton we find a good agreement up to 
the crossover time for the DST solution. Then the DST exciton can be seen to decay faster 
than the polaron (0). The reason is that the initially localized state of the exciton has for 
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the two solutions different sources. For the DST solution it results from self trapping on 
the attractive fixed point which keeps the exciton localized as long as it is far from its 
threshold of existence. Once this threshold is reached, the exciton becomes delocalized for 
good. 

In contrast, for the polaron solution the oscillator is too inert to change its position and 
thus keeps the exciton in a fixed adiabatic state. While B~ represents a unique point on 
the Bloch sphere for given parameters and total occupation, there is nothing special about 
the adiabatic self trapped state. In fact the exciton could for any set of parameters be fixed 
anywhere if the initial condition for the oscillator was chosen appropriately. Moreover, as 
the oscillator changes slowly its position, there will always be intervals when the exciton is 
localized on the sink site, i. e. the initial localization of the exciton is not due to self trapping 
and the initial agreement of the DST solution and the polaron (0) is simply due to the fact 
that the system was prepared in a DST state. 

This interpretation is further confirmed by the observation that the initial agreement 
between the DST and the polaron (0) solution ceases to exist as soon as the parameters do 
not support a self trapped state for the DST case. The polaronic solutions are unaffected 
by this and do still display an initial tendency towards localization on the sink less site (not 
displayed). Among them the polaron (+) again decays fastest while the polaron (-) is the 
longest living solution. 

The situation is similar in fig. |^. The DST exciton relaxes due to the strong sink quite 
fast to the final location of the fixed point B~ and then decays without further oscillations 
while the polaron (0) due to its inertness remains for a longer time close to its initial position 
and decays consequently slower than the DST solution. All the polarons as well as the bare 
exciton keep oscillating around a mean value which is slightly below the location of the 
fixed point B~. In fact fig. |^(b) looks very much like the plots for (7 = 3 and p = 3 at 
high and intermediate frequency fig. 0(b) and fig. ^(b), just the deviation from the location 
of the point B~ is larger and the oscillations are slower. But now we can provide a more 
satisfactory explanation for this behavior using the adiabatic states. The polarons as well 
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as the bare exciton follow after a very short relaxation the adiabatic state at cos 6 < 0. As 
an example for this behavior the adiabatic state for the polaron (-) is displayed in fig. ^(b) 
with sparse fat dots. The location of the adiabatic state can be seen from ( |39| ) to depend on 
the squared amplitude of the oscillator coordinate. The localization is weakest for Q = 0, 
when the adiabatic state coincides with the point B^. The quadratic rather than linear 
dependence on Q is the reason why the adiabatic oscillations in the relative site occupation 
observed in fig. ^(b) have a mean value below and a frequency which is exactly half that 
of the oscillator. 

In contrast to all the other solutions, the bare exciton remains completely unaffected 
by the nonlinearity in the adiabatic case. Here, the oscillator is prepared at Q = and 
there it stays during the whole life time of the excitation provided the adiabatic parameter 
r is sufficiently small. Consequently the vibronic coupling has no effect on the exciton and 
it oscillates independent on p around cos 6 = for < 1 (fig. or relaxes to cos 6' = 
^1 — otherwise (fig. |^). When the coupling parameter p is small, the bare exciton is 
well approximated by the DST solution. In the adiabatic regime r <^ 1 the bare exciton 
shows among the different considered solutions at least initially the fastest decay. 

Finally we would like to discuss an example in which the condition ( |5(]| ) for an exponential 
decay of the polaron solutions is satisfied (fig. §). In this case the polarons do not differ 
very much from each other and clearly follow an exponential law at a rate very close to 
that predicted by (|49|). Since we chose p = g for the figure, the life time of the polarons is 
exactly twice that of the bare exciton. The little remaining difference between the polaronic 
solutions reflects the residual change in the oscillator position during the life time of the 
excitation which enhances the localization of the exciton for the polaron (-) and diminishes 
it for the polaron (+) while there is no such effect for the polaron (0). 
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V. CONCLUSIONS 



We have studied the decay of an exciton coupled to polarization vibrations on a dimer. 
Quasistationary decay modes were identified which allow to explain the basic properties of 
the system. Using numerical simulations the deviations from the predicted behavior were 
investigated. 

The model exhibits a rich variety of dynamical regimes depending on the parameters 
and the initial conditions. We found effects such as the time dependent bifurcation and the 
associated crossover in the decay regime which are genuinely due to the interplay between 
the sink and the vibrational coupling and cannot be explained by considering one of these 
mechanisms alone. 

The tendency to form an initially localized exciton state on the site without sink is 
enhanced by both, vibrational coupling and trapping due to the sink. For high and inter- 
mediate oscillator frequency the system changes its behavior profoundly when the threshold 
for an initially self trapped state is reached, while there is no such effect in the adiabatic 
regime. 

The relation between the DST approximation and a mixed quantum-classical description, 
taking the oscillator dynamics explicitly into account, was clarified. For high oscillator 
frequency the infiuence of the oscillator initial condition is weak and the two models behave 
very much the same. In the adiabatic regime the bare exciton is close to the DST solution 
provided that the coupling is weak. 

The strong dependence on the initial conditions in the adiabatic case make a careful 
description of the exciton creation process indispensable for a satisfactory description of the 
system. Other interesting possibilities to extend the model are the inclusion of dissipation 
and / or quantum fiuctuations. 
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FIG. 1. Time dependence of the total occupation (top) and of the relative site occupation 
difference (bottom) for 5 = 0.1, p = 1 and r = 10. Different oscillator initial conditions for the full 
dynamic model and the DST dynamics can hardly be distinguished for this parameter set. They 
arc shown in this and all the following figures with the Unc types indicated in the upper part. In 
the bottom plot, cos^ = 1 corresponds to the sink site and cos^ = — 1 to the sinkless site where 
the exciton is created. The lower/upper horizontal line shows the location of the fixed point B"~ at 
the time of the creation of the exciton {R = 1) and after its complete decay {R = 0), respectively. 
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FIG. 3. Time dependence of the total occupation (top) and of the relative site occupation 



difference (bottom - DST, bare exciton and polaron (0) only) for = 3, p = 3 and r = 10. 



36 




37 




38 



0.01 



p' 

o 
o 



10 20 30 40 50 60 70 80 90 100 

X 



1 

0.8 
0.6 
0.4 

0.2 


-0.2 
-0.4 
-0.6 
-0.8 
-1 



iyiJiJVfe.--^;;........!.. 



10 20 30 40 50 60 70 80 90 100 

X 

FIG. 6. Time dependence of the total occupation for all initial conditions (top) and of the 
relative site occupation difference for the polaron (-) (bottom). With the sparse bold dots the 
time dependence of the two adiabatic states is indicated. The parameters are g = 0.1, p = 3 and 
r = 0.1. 
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FIG. 7. Time dependence of the total occupation (top) and of the relative site occupation 
difference (bottom) ioi g = 3, p = 3 and r = 0.1. In the bottom plot for the polaron (-) solution 
beside the relative site occupation the time dependence of the lower adiabatic state is displayed 
with sparse bold dots. 
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FIG. 8. Time dependence of the total occupation for g = 3, p = 3 and r = 0.001. 
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